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Abstract 

The dispersive character of the Hall-MHD solutions, in particular the whistler waves, 
is a strong restriction to numerical treatments of this system. Numerical stability 
demands a time step dependence of the form At oc (Ax)^ for explicit calculations. A 
new semi-implicit scheme for integrating the induction equation is proposed and ap- 
plied to a reconnection problem. It it based on a fix point iteration with a physically 
motivated preconditioning. Due to its convergence properties, short wavelengths 
converge faster than long ones, thus it can be used as a smoother in a nonlinear 
multigrid method. 
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1 Introduction 



In many space-, astrophysical and high temperature plasma systems collisions 
do not play the most important role in describing the departure from the ideal 
magnetohydrodynamics (MHD) 

dtp = -V • [pv) (1) 

dtv = -{v-W)v + ^^-^ (2) 

P P 

dtB = -V X ^ (3) 

j = VxB, (4) 

where p, v, B and p denote mass density, velocity, magnetic field and pres- 
sure, respectively. Typical examples include filamentation and singularity for- 
mation, coUisionless reconnection and coUisionless shocks [1,2,3,4,5,6,7,8,9]. 
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Therefore, on scales smaller than the ion inertia length additional processes 
have to be taken into account in a generalized Ohm's law 



E^rjj-vxB+^(jxB-Vpe) . (5) 
Numerically, the most difficult term is the Hall-term: 

EhoXI = -w^J xB = ^jxB. 6 
Zep p 

It allows for whistler wave solutions with a quadratic dispersion relation and 
thus poses a severe time step restriction for a temporal explicit discretisation. 

To introduce our treatment of the Hall-term, we simplify our system and use 
only this electric field in the induction equation which decouples it from the 
other part of the MHD equations and yields the following nonlinear equation 



dtB = -V X (V X 5) X 5 j . 



(7) 



Solutions of the linearized equations are the whistler waves mentioned above 
which satisfy the dispersion relation u — ^^fc^, for a constant density p and 

a guiding field magnitude \B\. Numerical approaches using explicit schemes 
applied to this equation must ensure that the chosen time step fulfills At oc 
{Axy, due to the Courant-Levy-Friedrichs criterion - Ax denoting the grid 
spacing. The CFL number is given by the ratio of the phase velocity to the 
grid velocity (^^^ 



B 



k Ax pAx^ di B 

where k — k^ax — ^ maximum wave number. Thus resolving small 

structures, e.g. the reconnection zone, results in large computation times, due 
to the unavoidable small time steps. 

Implicit schemes allow to avoid this restrictive condition by providing uncon- 
ditional numerical stability. Much progress on implicit solvers has been done 
by Harned and Mikic [10] and Chacon and Knoll [11]. However, the approach 
of [10] requires a guiding magnetic field and the approach of [11] can't eas- 
ily be adopted for simulations with adaptive mesh refinements [12,13,14,15], 
although work in this direction is in progress. 

Here we present a simple physics based semi-implicit Crank-Nicolson type 
scheme which due to its locality properties is suitable for parallel computations 
as well as for use in adaptive mesh refinement simulations. This physics based 



2 



solver Tiscs a whistler wave decomposition to accelerate the fix-point iteration. 
Due to its convergence properties it can act as a smoother for a nonlinear 
multigrid scheme. 

The first part of this paper presents the general numerical method which then 
is specialized to one dimension. This allows us to show analytically its con- 
vergence. After that the nonlinear two-dimensional case and its convergence 
are presented, while in the last section our method is used to solve a two- 
dimensional reconnection problem. 



2 Numerical Method 



The Richardson iteration [16] is the base of our solver. A Richardson iteration 
is the most general fix point iteration for a nonlinear equation F{x) — 

x''+^^K{x'') with K{x) = x-aF{x) , (8) 

where k is the iteration index. Given a contractive map K, the x'^ converge in 
the limit k oo. The rate of convergence will in general depend on a. The 
main task is to find a suitable preconditioner adapted to the Hall-term. This 
can be realized as a matrix P 

K{x) = f - aPF{x) . (9) 

In the special case of the Newton iteration P is the inverse of the Jacobi matrix 
of F. Here, we try to find a physics based preconditioner which is more local 
and thus suitable for parallel and block-adaptive calculations. 

For the Crank- Nicolson type discretisation, we obtain 



At 




with B* = 1(5"+^ -I- S") and where is the magnetic field taken at the time 
step n (time steps are indicated by the first upper index) . The equation to be 
solved reads now 

on+l on / J \ 

^ +Vxi-(yxB*)xB*\ (11) 

= 0. (12) 

Its solution with a given fi" is the magnetic field at the next time step n + 1. 
To determine a solution we iterate equation (11) following the method given 
by (8). At this point we introduce an additional upper index which defines 
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the iteration step. So that is k-ih iteration of the magnetic field for the 

time step n + 1. 

As mentioned above, the important point in this iteration is the precondition- 
ing. To motivate our preconditioner, we start with the one-dimensional version 
of eqn. (7), where B depends only on x. In the one-dimensional case, we have 
the special situation that the fix-point problem reduces to a linear one. In 
this case, our physics based preconditioner reduces to the standard cj-Jacobi 
iteration. In more than one dimensions the situation is genuinely nonlinear 
but the smoothing properties of our preconditioner are still similar to that of 
a Jacobi iteration for linear problems. 



3 Linear ID Ccise 



Considering only the x direction results in the following system of equations 



By 

yB^j 



\ 



^dxx{B,) 

-^dxABy) ) 



(13) 



The X component of B is initially set to a constant Bq in space and stays con- 
stant. The discretised equations for F arc needed for the iteration. Following 
equation (11) and again using B* = ^ (^B"' + these are given for each 

grid point i 



F.{Br') =5^' - b: 

?n+l 



X,l 



Fy{Br') ^B^r - - c.{b:,_, - 2b:, + 5;,^,) 

F,(Br') ^b:^' - Bi, + c.(5^_, - 2b;. + b;.^,) 



(14) 



with Cx = ^^^0- To calculate the next time step B'^'^^ we consider the itera- 
tion such that hmfe^oo 5"+^''= = 



The usual choice for a preconditioning matrix would be the inverse of the 
full Jacobian matrix. Even in the one-dimensional case this would lead to 
an inversion of a 2N x 2N matrix (A^ is the number of grid points). In our 
treatment of the Hall-term, we introduce a local approximation to the Jacobian 
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matrix. The Richardson iteration (8) for F = reads 





=Bx,i 












■r)n+l,k 














+1' 




— 2B^~^^'^ 


1 r>n+l,k 
+ ^y,i+i 



(15) 



where rx,i,ry^i,rz,i are constants depending only on the values B"'. 

The matrix elements of the local preconditioner are calculated by differenti- 
ating Eqs. (15) with respect to B^'^^''^; this is a differentiation with respect to 
the value of B at one grid point, e.g. the derivative of By^^'i with respect to 
By^^''' vanishes. This leads to the following local preconditioner 



'^l ^ 

1 Cr 



-c^ 1 



(16) 



The preconditioning matrix J* is now a block diagonal matrix containing only 
the J*. This local construction allows an easy inversion, which is again a block 
matrix and thus (J*)~^ is local. A fast inversion of the general Jacobian is 
not easy and results in a non local matrix. Using (J*)~^ as a preconditioner 
results in the following local iteration for each grid point i. 



^n+l,k+l 



a 



l + cl 



+ ^ 
1 -c, 
1 



(17) 



This iteration is the cu- Jacobi iteration, which in this case is known to converge. 
This way we have created a well known iteration scheme, but it was motivated 
by the physical properties of a whistler wave. While the Jacobi iteration cannot 
be used for nonlinear problems, we can transfer our iteration to the nonlinear 
two-dimensional case using the same strategy. The iteration obtained so far 
has the desirable property that high wavenumbers converge fast and thus can 
be used as a smoother in a multigrid scheme. We anticipate that this property 
translates to the two-dimensional case. 
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4 Nonlinear 2D case 



The strategy here will be the same as above: 



i,3 hi 



(18) 



The local preconditioner is motivated by the one-dimensional calculation, tak- 
ing into account the two directions x and y of whistler wave propagation. In 
two dimensions the function F can be derived the same way as the in the 
one-dimensional case and takes the form 



dyB;dyB; + d^Btd^Bt + Bld^^B 



Fy = - 5; - 



B 



n+l 



b: + 



p 

diAt 

P 

diAt 



{B;d,yB: + d,B;dyB: + d^Bld^Bl + Bld^^Bl) 



{d^Bi + {d^B; - dyBi) (19) 

+ -B* {d:cxB*y - d^yB*}j + B* {dxyB*y - dyyBf^ 



To obtain the local preconditioner, Eq. (19) has to be discretised in space and 
derived with respect to B^^^ which yields 



/ 



T* — 



1 









1 


Cx,i,j 


Cx, 


ij 1 



(20) 



'-^ — and Cy^ij 



— _ Note that the elements of this 



with C^^ij = pA:r^ ^^^^ = ~7aP 

Jacobian matrix are not constant anymore, but depend on -B^^jj and B^ j^j 
Thus again the preconditioning is it's inverse 



"I" ^x,i,j ^x,i,jCy,i,j Cy^ij 



T*— 1 



^ y,i,j 



Cx,i,jCy,i,j 



V 



(21) 



To show numerically the convergence properties of this iteration, a simple 
numerical experiment is used. The two-dimensional computational domain is 
a periodic box in x- and y-direction. The initial condition is a single whistler 
wave with a wave number in x- and y-direction. This is done for all possible 
modes (combinations of and ky ) on a 64 x 64 mesh. The number of iterations 
n eeded fo r a prescribed accuracy is plotted in Fig. 1 as a function of A; = 
Jkl + k^. As we already anticipated the convergence rate for the long and the 
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Fig. 1. The number of iterations needed to reduce the error by two orders of mag- 
nitude as a function of the norm of the wave vector k = ^k"^ + ky. 

short wavelengths show the same behavior as in the one-dimensional case. 

To accelerate the convergence of the long wavelengths, this iteration is applied 
as a smoothing function for the nonlinear multigrid scheme [17]. Using V- 
cycles and two pre- and post-smoothings in the multigrid scheme one achieves 
convergence already after two to three cycles. 



5 GEM Reconnection 



To verify the new iteration scheme, we choose a standard reconnection prob- 
lem. A similar setup is used as in GEM reconnection challenge [7]. It is based 
on a perturbed Harris sheet. A schematic plot and the computational domain 
are shown in Fig. 2. The numerical parameter are chosen to: = 2Ly = 5di 
and Nx = 2Ny = 256. Symmetric or antisymmetric boundary conditions are 
applied for all quantities and the initial conditions, including the perturbation, 
are equal to the one in the GEM reconnection challenge [7]. 

We solve the full Hall-MHD equations explicitly except for the Hall part of 
the induction equation. This splits the induction equation into the resistive 
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Fig. 2. The GEM reconnection setup. The black hnes indicate the magnetic field, 
the blue ones the seperatrices and the red box is our computational domain. 



MHD part (22) and the part including only the Hall term (23) 





(22) 



(23) 



dtB 



(24) 



The time stepping was performed with a standard second order Runge-Kutta 
method. 

The maximum time step for the semi-implicit simulations is given by the 
linear Alfven wave dispersion relation uo = vjjt. There is no possibility to use 
larger time steps than these, because the Alfven waves must be well resolved. 

To compare the results obtained with the fully explicit simulation and with our 
semi-implicit treatment of the Hall term, we choose as a physically relevant 
measure the reconnected flux ip given in our setup by 



Four different simulation runs have been done, an explicit one, used as the 
reference run, and three semi-implicit ones. The time steps and corresponding 
CFL numbers are shown in Table 1. 




(25) 



Fig. 3 shows the developed structure of the z-component of the electric current 
density j = V x i? at time t = 12 obtained from a semi-implicit simulation. 





At ■ 10"^ 


riT^'Ti TniTTiHpr 


^^exp 




0.2 


0.2 


1 


semi-implicit 


2.0 


2 


10 


semi-implicit 


4.0 


4 


20 


semi— implicit 


8.2 


8.2 


41 



Table 1 

Time steps chosen for the simulations. The CFL number is based on the whistler 
wave dispersion relation. 

The corresponding reconnected flux (Fig. 4) obtained from the explicit refer- 




Fig. 3. The z-component of the electric current density at t = 12. 

ence simulation and the semi-implicit runs show that the the semi-implicit 
simulations result in nearly identical reconnection rates and that, by reducing 
the time step, they converge to the values obtained from the explicit simula- 
tion. 



6 Summary 

Wc presented a semi-implicit iterative method for solving the Hall part of 
the induction equation to overcome the time step restriction resulting from 
the quadratic whistler wave dispersion relation. The method utilizes a simple 
precondioner based on the whistler wave dispersion. The iteration scheme 
based on this precondioner has the two desirable features of being local and 
possessing strong high frequency smoothing. Therefore, this method can easily 
be implemented in a nonlinear multigrid solver. Due to the locality property, it 
is also best suited for adaptive and parallel simulations. The part of execution 
time of Hall term turns out to be about 7 times longer than the time needed 
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Fig. 4. The temporal evolution of the reconnected flux for different time steps and 
schemes. 

for ideal MHD part. Using a time step 40 times larger than necessary for an 
explicit treatment, this results in an 80% reduction of the computation time 
for the GEM setup achieving nearly identical results as an expensive explicit 
simulation. 
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